!!!! Business Cycles, Insurance - Deterministic SS Model - (3/11/2024) !!!!
!! Author: Michael Irwin
!! from "The Impact of Unemployment Insurance and Unsecured Credit on Business Cycles"
!! for JPE Macroeconomics
!! Solution Method (individuals): OEGM for savings; grid search for debt
!! Solution Method (aggregate): Steady-State Equilibrium

!!!! State Variables !!!!
!! e : Idiosyncratic Productivity
!! a : Asset Level
!! j : Age
!! s : Credit Status
!! n : Individual Labor Market Status
!! b : Beta (discount factor)
    
!!!! Code Options !!!!
!! 1.) 'switch_scenario', determines which scenario is run
!! 2.) 'switch_outerloop', solves for a fixed point for K:L ratio
!! 3.) 'switch_welfare' collects Vc to calculate welfare
    
!!!! Scenarios !!!!
!! 1.) Benchmark: parameters from 'expansion' state in stochastic equilibrium model
!! 2.) Recession analysis: Parameters for recession
!! 3.) Recession analysis: recession value for TFP
!! 4.) Recession analysis: recession value for job separation rates
!! 5.) Recession analysis: recession value for job finding rates
!! 6.) UI Extension: extend duration of UI
!! 7.) UI RR Increase: increase RR of UI
!! 8.) No Credit: model with no unsecured credit
!! 9.) No Credit, UI Extension
!! 10.) No Credit, UI RR Increase
!! 11.) Recession, UI extension
!! 12.) Recession, RR Increase
!! 13.) Recession, No Credit
!! 14.) Recession, No Credit, UI extension
!! 15.) Recession, No Credit, RR Increase
!! 16.) No UI: Model with replacement rate set to 0.0
!! 17.) No Credit, No UI

!! Code Introduction
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Model Parameters

program BCI_SS


use BCI_Par
use BCI_TB
use OEGM_TB
use Fortran_TB
use omp_lib

implicit none

!!!! Switches !!!!
integer, parameter :: switch_scenario=14, switch_outerloop=1, switch_results=1, switch_welfare=1

!!!! Variables for Solving Model !!!!
integer :: iteration, count, seedsize
integer :: a, e, j, n, s, b, i, p, t, a2, e2, n2, i2, ia, ij, ibar_e, ibar_u, ip
integer :: ih, ie, ix, it
integer :: ngrids_ge, ngrids_gu, ngrids_be, ngrids_bu
real(rp) :: error, omega
real(rp) :: rate, wage, rate_save
real(rp) :: Value, V_d, y_u, chi_e, chi_u, income_e, income_u
real(rp) :: maxV_e, maxG_e, maxC_e, maxD_e, maxQ_e, maxV_u, maxG_u, maxC_u, maxD_u, maxQ_u
real(rp) :: muval, gval, dval, hval, cval
real(rp) :: avg_earn
real(rp) :: abar_e, abar_u
real(rp) :: lambda, xi, zgrid, upsilon_d, upsilon_r, upsilon_bar, iota
real(rp) :: r, w, E_r, aggK, aggKp, aggL, aggLp, error_K, error_L

!!!! Variables for Results !!!!
integer, parameter :: n_income=nj*ni
integer, dimension(n_income) :: i_income, j_income
real(rp), dimension(n_income) :: income_grid, mu_income
real(rp) :: aggY, aggC, aggI, aggG, aggD, aggB, aggXi
real(rp) :: Dpop, Upop, Epop, Npop
real(rp) :: aggC_l, aggC_h, aggC_debt, aggC_save
real(rp) :: Dpop_E, Dpop_U, Dpop_P, Dpop_N
real(rp) :: UI_exp, ret_exp, T_exp, tax_rev
real(rp) :: sum_mu, qval, aggA, agg_discharge, profits

!!!! Variables for Borrowing Limits and Debt Prices !!!!
real(rp) :: avg_q, avg_limit, pop_borrow, blim, limit_all, limit_patient, limit_impatient, pop_all, pop_patient, pop_impatient

!!!! Arrays !!!!
integer, dimension(ni) :: n_index, e_index, b_index
integer, dimension(no) :: ilow_ge, ihigh_ge, ilow_gu, ihigh_gu, ilow_be, ihigh_be, ilow_bu, ihigh_bu
integer, dimension(ne,nn,nb) :: i_index
real(rp), dimension(na) :: agrid, EnGrid_ge, EnGrid_gu, EnGrid_be, EnGrid_bu
real(rp), dimension(na) :: EV_ge, EV_gu, EV_be, EV_bu, EVc_ge, EVc_gu, EVc_be, EVc_bu, DV_ge, DV_gu, DV_be, DV_bu
real(rp), dimension(nd) :: tempgrid, discount_e, discount_u
real(rp), dimension(ne) :: egrid, pi_erg, cum_pi_erg
real(rp), dimension(nj) :: gamma_j
real(rp), dimension(nb) :: beta, phi_b
real(rp), dimension(ne,ne) :: pi_e, cum_pi_e
real(rp), dimension(nn,nn) :: pi_n, cum_n
real(rp), dimension(ni,ni,2) :: pi_i
real(rp), allocatable :: mu(:,:,:,:), mu_next(:,:,:,:), mu_good(:,:,:,:), mu_bad(:,:,:,:), mur(:,:,:), mur_next(:,:,:)
real(rp), allocatable :: V(:,:,:,:), V_c(:,:,:,:), g(:,:,:,:), c(:,:,:,:), d(:,:,:,:), h(:,:,:,:), Vr(:,:,:), gr(:,:,:), cr(:,:,:), abar(:,:,:), q(:,:,:)
integer, allocatable :: seed(:)

!!!! Set the Number of Threads !!!!
!$ call omp_set_num_threads(np)

!!!! Variables for Aggregate Fluctuations !!!!
iota = iota_g
if ( switch_scenario == 1 .or. switch_scenario == 6 .or. switch_scenario == 7 .or. switch_scenario == 8 .or. switch_scenario ==9 .or. switch_scenario == 10 .or. switch_scenario == 16 .or. switch_scenario == 17) then
    zgrid = z_g; lambda = lambda_g; xi = xi_g
elseif ( switch_scenario == 2 .or. switch_scenario == 11 .or. switch_scenario == 12 .or. switch_scenario == 13 .or. switch_scenario == 14 .or. switch_scenario == 15 ) then
    zgrid = z_b; lambda = lambda_b; xi = xi_b
elseif ( switch_scenario == 3 ) then
    zgrid = z_b; lambda = lambda_g; xi = xi_g
elseif ( switch_scenario == 4 ) then
    zgrid = z_g; lambda = lambda_b; xi = xi_g
elseif ( switch_scenario == 5 ) then
    zgrid = z_g; lambda = lambda_g; xi = xi_b
end if

!!!! Variables for UI Policies !!!!
upsilon_bar = upsilon_bar_g
if ( switch_scenario <= 5 .or. switch_scenario == 8 .or. switch_scenario == 13 ) then
    upsilon_d = upsilon_dg; upsilon_r = upsilon_rg
elseif ( switch_scenario == 6 .or. switch_scenario == 9 .or. switch_scenario == 11 .or. switch_scenario == 14) then			!! UI Extension
    upsilon_d = upsilon_db; upsilon_r = upsilon_rg
elseif ( switch_scenario == 7 .or. switch_scenario == 10 .or. switch_scenario == 12 .or. switch_scenario == 15) then			!! UI RR Increase
    upsilon_d = upsilon_dg; upsilon_r = 0.62_rp
elseif ( switch_scenario == 16 .or. switch_scenario == 17 ) then        !! No UI (RR=0)
    upsilon_d = upsilon_db; upsilon_r = zero
end if

!!!! Grid of Financial Wealth !!!!
call logspace(zero, -a_min, nd, tempgrid)
do a = 1,nd
	agrid(a) = -tempgrid(nd+1-a)
end do
call logspace(zero, a_max, na-nd+1, agrid(nd:na))

!!!! Persistent Earnings Process !!!!
call rouwenhorst(rho, sigma_eta_g, ne, pi_e(:,:), egrid)
call ergodic(ne, pi_e(:,:), pi_erg(:))
egrid = exp(egrid)
cum_pi_e(:,1) = pi_e(:,1)
cum_pi_erg(1) = pi_erg(1)
do e = 2,ne
    cum_pi_e(:,e) = cum_pi_e(:,e-1) + pi_e(:,e)
    cum_pi_erg(e) = cum_pi_erg(e-1) + pi_erg(e)
end do

!!!! Age-Component of Earnings !!!!
do j = 1,nj,4
	ij = ceiling(j/4.0_rp)
	gamma_j(j:j+3) = nu_1*dble(ij) + nu_2*dble(ij)**2.0_rp
end do
gamma_j = exp(gamma_j)

!!!! Transition of Employment States !!!!
call emp_transitions(nn, xi, lambda, upsilon_d, pi_n(:,:))
cum_n(:,1) = pi_n(:,1)
do n = 2,nn
    cum_n(:,n) = cum_n(:,n-1) + pi_n(:,n)
end do

!!!! Heterogeneous Discount Factors !!!!
beta(1) = beta_h
beta(2) = beta_l
phi_b(1) = pi_beta
phi_b(2) = one-pi_beta

!!!! Assign Indices and Create Transition Matrix !!!!
call ind_transitions(ni, ne, nn, nb, pi_e, pi_n, e_index, n_index, b_index, i_index, pi_i(:,:,:) )

!!!! Set Income Distribution !!!!
call income_distribution(n_income, e_index, n_index, upsilon_r, gamma_j, egrid, i_income, j_income, income_grid)

!!!! Allocate Large Arrays !!!!
allocate(V(ni,na,nj,ns), g(ni,na,nj,ns), c(ni,na,nj,ns), d(ni,na,nj,ns), h(ni,na,nj,ns))
allocate(Vr(nb,na,nR), gr(nb,na,nR), cr(nb,na,nR))
allocate(abar(ni,na,nj), q(ni,na,nj))
allocate(mu(ni,na,nj,ns))
allocate(mur(nb,na,nR))
if ( switch_welfare == 1 ) then
    allocate(V_c(ni,na,nj,ns))
end if

!!Parameters of the Model
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!Begin Outer Loop

!!!! Outer Loop - Convergence of K:L Ratio !!!!
error = 100.0_rp
iteration = 0
aggK = 38.28_rp
aggL = 1.31_rp
!aggKL = 32.0_rp

do while ( error >= 0.0001_rp .and. iteration < 500 )

    !!!! Prices !!!!
    rate = alpha*zgrid*( aggK**(alpha-one) )*( aggL**(one-alpha) ) - delta
    wage = (one-alpha)*zgrid*( aggK**alpha )*( aggL**(-alpha) )
    !rate = alpha*zgrid*( aggKL**(alpha-one) ) - delta
    !wage = (one-alpha)*zgrid*( aggKL**alpha ) 
    rate_save = rate

!!Begin Outer Loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!Solving for Decision Rules 

!!!! Initial Functional Values !!!!
V = zero; g=zero; c=zero; d=zero; h=zero; Vr=zero; gr=zero; cr=zero
abar=zero; q=zero
if ( switch_welfare == 1 ) then
    V_c = zero
end if

!!!! Final Age Value !!!!
do a = nd,na
    Vr(:,a,nR) = (one/(one-sigma))*(( agrid(a) + y_r )**(one-sigma))
    cr(:,a,nR) = agrid(a) + y_r
    gr(:,a,nR) = zero
end do

!!!! Decision Rules in Retirement - EGM !!!!
do j = (nR-1),1,-1
	do b = 1,nb

        !!!! Derivatives of Expected Future Value !!!!
    	EV_gu = zero
		do a2 = na,nd,-1
            EV_gu(a2) = Vr(b,a2,j+1)

            if ( a2 < na ) then
        		DV_gu(a2) = ( EV_gu(a2+1)-EV_gu(a2) )/( agrid(a2+1)-agrid(a2) )
            end if
		end do
        DV_gu(na) = DV_gu(na-1)

        !!!! Solve for Decision Rules !!!!
        call endogenous_grids(beta(b), rate_save, y_r - y_bar, agrid, DV_gu, ngrids_gu, ilow_gu, ihigh_gu, EnGrid_gu)
        maxV_u = disutility; maxG_u = zero; maxC_u = zero; maxQ_u = zero
        do a = nd,na
            call oegm_interpolation(ngrids_gu, beta(b), zero, agrid(a), y_r - y_bar, rate_save, ilow_gu, ihigh_gu, EnGrid_gu, agrid, EV_gu, maxV_u, maxG_u, maxC_u, maxQ_u)
            Vr(b,a,j) = maxV_u
        	gr(b,a,j) = maxG_u
    		cr(b,a,j) = maxC_u
        end do

	end do
end do

!!!! Decision Rules Working Ages - OEGM !!!!
do j = nj,1,-1

	!$omp parallel do private(i, a, a2, n, e, b, s, discount_e, discount_u, ibar_e, ibar_u, abar_e, abar_u, &
    !$              EV_ge, EV_gu, EV_be, EV_bu, EVc_ge, EVc_gu, EVc_be, EVc_bu, DV_ge, DV_gu, DV_be, DV_bu, EnGrid_ge, EnGrid_gu, EnGrid_be, EnGrid_bu, &
    !$              ilow_ge, ilow_gu, ilow_be, ilow_bu, ihigh_ge, ihigh_gu, ihigh_be, ihigh_bu, ngrids_ge, ngrids_gu, ngrids_be, ngrids_bu, &
	!$              maxV_e, maxG_e, maxC_e, maxD_e, maxQ_e, maxV_u, maxG_u, maxC_u, maxD_u, maxQ_u, p, ia, omega, &
	!$              Value, V_d, y_u, income_e, income_u, chi_e, chi_u )
	

	do i = 1,ni
		n = n_index(i); e = e_index(i); b = b_index(i) 
            
            !!!! Variables Dependent on Labor Choice !!!!
            if ( n == 1 ) then
                income_e = (one-tau)*wage*egrid(e)*gamma_j(j)
                income_u = (one-tau)*wage*egrid(e)*gamma_j(j)
                chi_e = chi_w
                chi_u = chi_w
            elseif ( n == 2 ) then
                if ( upsilon_r*wage*egrid(e)*gamma_j(j) >= upsilon_bar ) then
                    income_e = upsilon_bar
                else
            		income_e = upsilon_r*wage*egrid(e)*gamma_j(j)
                end if
                income_u = zero
                chi_e = chi_s
                chi_u = zero
            else
                income_e = zero
                income_u = zero
                chi_e = chi_s
                chi_u = zero
            end if
            
            !!!! Calculate Expected Future Variables !!!!
            if ( j == nj ) then
                
                discount_e(:) = zero; discount_u(:) = zero
                ibar_e = nd; ibar_u = nd
                EV_ge = zero; EV_gu = zero; EV_be = zero; EV_bu = zero
                if ( switch_welfare == 1 ) then
                    EVc_ge = zero; EVc_gu = zero; EVc_be = zero; EVc_bu = zero
                end if
                
		        do a2 = na,nd,-1
                    EV_ge(a2) = Vr(b,a2,1); EV_gu(a2) = Vr(b,a2,1); EV_be(a2) = Vr(b,a2,1); EV_bu(a2) = Vr(b,a2,1)
                    if ( switch_welfare == 1 ) then
                        EVc_ge(a2) = Vr(b,a2,1); EVc_gu(a2) = Vr(b,a2,1); EVc_be(a2) = Vr(b,a2,1); EVc_bu(a2) = Vr(b,a2,1)
                    end if

                    if ( a2 < na ) then
                	    DV_ge(a2) = ( EV_ge(a2+1)-EV_ge(a2) )/( agrid(a2+1)-agrid(a2) )
            	        DV_gu(a2) = ( EV_gu(a2+1)-EV_gu(a2) )/( agrid(a2+1)-agrid(a2) )
            	        DV_be(a2) = ( EV_be(a2+1)-EV_be(a2) )/( agrid(a2+1)-agrid(a2) )
                	    DV_bu(a2) = ( EV_bu(a2+1)-EV_bu(a2) )/( agrid(a2+1)-agrid(a2) )
                    end if
                    
    	        end do
                DV_ge(na) = DV_ge(na-1); DV_gu(na) = DV_gu(na-1); DV_be(na) = DV_be(na-1); DV_bu(na) = DV_bu(na-1)
                
            else
                
                !!!! Discount Prices !!!!
			    discount_e(:) = (one-iota)/(one+rate_save)
                discount_u(:) = discount_e(:)
                do a2 = 1,(nd-1)
                    discount_e(a2) = discount_e(a2) - dot_product(pi_i(i,:,1),d(:,a2,j+1,1))/(one+rate_save)
                    discount_u(a2) = discount_u(a2) - dot_product(pi_i(i,:,2),d(:,a2,j+1,1))/(one+rate_save)
                end do

			    !!!! Calculate the Endogenous Borrowing Limit !!!!
                if ( switch_scenario == 8 .or. switch_scenario == 9 .or. switch_scenario == 10 .or. switch_scenario == 13 .or. switch_scenario == 14 .or. switch_scenario == 15 .or. switch_scenario == 17 ) then
                    ibar_e = nd
                    ibar_u = nd
                else
			        call borrowing_limit(nd, agrid(1:nd), discount_e(:), ibar_e, abar_e)
			        call borrowing_limit(nd, agrid(1:nd), discount_u(:), ibar_u, abar_u)
                end if
            
                !!!! Calculate the Expected Future Value and Numerical Derivatives !!!!
                EV_ge = zero; EV_gu = zero; EV_be = zero; EV_bu = zero
                do a2 = na,1,-1	
                    EV_ge(a2) = dot_product(pi_i(i,:,1),V(:,a2,j+1,1))
                    EV_gu(a2) = dot_product(pi_i(i,:,2),V(:,a2,j+1,1))
                    EV_be(a2) = dot_product(pi_i(i,:,1),V(:,a2,j+1,2))
                    EV_bu(a2) = dot_product(pi_i(i,:,2),V(:,a2,j+1,2))

        		    if ( a2 < na ) then
        			    DV_ge(a2) = ( EV_ge(a2+1)-EV_ge(a2) )/( agrid(a2+1)-agrid(a2) )
        			    DV_gu(a2) = ( EV_gu(a2+1)-EV_gu(a2) )/( agrid(a2+1)-agrid(a2) )
        			    DV_be(a2) = ( EV_be(a2+1)-EV_be(a2) )/( agrid(a2+1)-agrid(a2) )
        			    DV_bu(a2) = ( EV_bu(a2+1)-EV_bu(a2) )/( agrid(a2+1)-agrid(a2) )
        		    end if
        	    end do
        	    DV_ge(na) = DV_ge(na-1); DV_gu(na) = DV_gu(na-1); DV_be(na) = DV_be(na-1); DV_bu(na) = DV_bu(na-1)
            
                !!!! Expected Value of Consumption - Welfare Analysis !!!!
                if ( switch_welfare == 1 ) then
                    EVc_ge = zero; EVc_gu = zero; EVc_be = zero; EVc_bu = zero
                    do a2 = na,1,-1	
                        EVc_ge(a2) = dot_product(pi_i(i,:,1),V_c(:,a2,j+1,1))
                        EVc_gu(a2) = dot_product(pi_i(i,:,2),V_c(:,a2,j+1,1))
                        EVc_be(a2) = dot_product(pi_i(i,:,1),V_c(:,a2,j+1,2))
                        EVc_bu(a2) = dot_product(pi_i(i,:,2),V_c(:,a2,j+1,2))	
                    end do
                end if
            
            end if

            !!!! Form Endogenous Grids !!!!
            call endogenous_grids(beta(b), rate_save, income_e, agrid, DV_ge, ngrids_ge, ilow_ge, ihigh_ge, EnGrid_ge)
            call endogenous_grids(beta(b), rate_save, income_u, agrid, DV_gu, ngrids_gu, ilow_gu, ihigh_gu, EnGrid_gu)
            call endogenous_grids(beta(b), rate_save, income_e, agrid, theta*DV_ge(:) + (one-theta)*DV_be(:), ngrids_be, ilow_be, ihigh_be, EnGrid_be)
            call endogenous_grids(beta(b), rate_save, income_u, agrid, theta*DV_gu(:) + (one-theta)*DV_bu(:), ngrids_bu, ilow_bu, ihigh_bu, EnGrid_bu)

            !!!! Solve for Decision Rules !!!!
            do a = 1,na

                !!!! Good Credit - Work/Search !!!!
                s = 1
                call grid_search(ibar_e, beta(b), chi_e, agrid(a), income_e, discount_e, agrid, EV_ge, maxV_e, maxG_e, maxC_e, maxQ_e)
                call oegm_interpolation(ngrids_ge, beta(b), chi_e, agrid(a), income_e, rate_save, ilow_ge, ihigh_ge, EnGrid_ge, agrid, EV_ge, maxV_e, maxG_e, maxC_e, maxQ_e)
                V_d = (one/(one-sigma))*(( income_e + y_bar )**(one-sigma)) - chi_e - chi_b + beta(b)*EV_be(nd)

                !!!! Bankruptcy Decision - Work/Search !!!!
                if ( V_d > maxV_e ) then   
                	maxV_e = V_d
                	maxG_e = zero
                	maxC_e = income_e + y_bar
                	maxD_e = one
                    maxQ_e = zero
                else
                	maxD_e = zero
                end if

                !!!! Good Credit - Leisure !!!!
                call grid_search(ibar_u, beta(b), chi_u, agrid(a), income_u, discount_u, agrid, EV_gu, maxV_u, maxG_u, maxC_u, maxQ_u)
                call oegm_interpolation(ngrids_gu, beta(b), chi_u, agrid(a), income_u, rate_save, ilow_gu, ihigh_gu, EnGrid_gu, agrid, EV_gu, maxV_u, maxG_u, maxC_u, maxQ_u)
                V_d = (one/(one-sigma))*(( income_u + y_bar )**(one-sigma)) - chi_u - chi_b + beta(b)*EV_bu(nd)

                !!!! Bankruptcy Decision - Leisure !!!!
                if ( V_d > maxV_u ) then
                	maxV_u = V_d
                	maxG_u = zero
                	maxC_u = income_u + y_bar
                	maxD_u = one
                    maxQ_u = zero
                else
                    maxD_u = zero
                end if

                !!!! Labor Market Decision !!!!
            	if ( maxV_e >= maxV_u ) then
        			V(i,a,j,s) = maxV_e
    				g(i,a,j,s) = maxG_e
                    c(i,a,j,s) = maxC_e
                	d(i,a,j,s) = maxD_e
            		h(i,a,j,s) = one
                    q(i,a,j) = maxQ_e
                    !abar(i,a,j) = agrid(ibar_e)
                    abar(i,a,j) = abar_e
                    if ( switch_welfare == 1 ) then
                        call interpolation(agrid, maxG_e, 1, na, na, ia, omega)
                        V_c(i,a,j,s) = (one/(one-sigma))*(( maxC_e )**(one-sigma)) + beta(b)*( omega*EVc_ge(ia) + (one-omega)*EVc_ge(ia+1) )
                    end if
                        
				else
                    V(i,a,j,s) = maxV_u
                	g(i,a,j,s) = maxG_u
            		c(i,a,j,s) = maxC_u
        			d(i,a,j,s) = maxD_u
    				h(i,a,j,s) = zero
                    q(i,a,j) = maxQ_u
                    !abar(i,a,j) = agrid(ibar_u)
                    abar(i,a,j) = abar_u
                    if ( switch_welfare == 1 ) then
                        call interpolation(agrid, maxG_u, 1, na, na, ia, omega)
                        V_c(i,a,j,s) = (one/(one-sigma))*(( maxC_u )**(one-sigma)) + beta(b)*( omega*EVc_gu(ia) + (one-omega)*EVc_gu(ia+1) )
                    end if
                    
        		end if

                !!!! Bad Credit Status !!!!
                s = 2
                if ( a >= nd )  then

                    !!!! Decision Rules from Working !!!!
                	maxV_e = disutility; maxG_e = zero; maxC_e = zero; maxQ_e = zero
            		call oegm_interpolation(ngrids_be, beta(b), chi_e, agrid(a), income_e, rate_save, ilow_be, ihigh_be, EnGrid_be, agrid, theta*EV_ge(:)+(one-theta)*EV_be(:), maxV_e, maxG_e, maxC_e, maxQ_e)

                    !!!! Start Decision Rules from Leisure !!!!
                	maxV_u = disutility; maxG_u = zero; maxC_u = zero; maxQ_u = zero
            		call oegm_interpolation(ngrids_bu, beta(b), chi_u, agrid(a), income_u, rate_save, ilow_bu, ihigh_bu, EnGrid_bu, agrid, theta*EV_gu(:)+(one-theta)*EV_bu(:), maxV_u, maxG_u, maxC_u, maxQ_u)

                    !!!! Labor Supply Decision !!!!
                	if ( maxV_e >= maxV_u ) then
            			V(i,a,j,s) = maxV_e
        				g(i,a,j,s) = maxG_e
                    	c(i,a,j,s) = maxC_e
                    	d(i,a,j,s) = zero
                    	h(i,a,j,s) = one
                        if ( switch_welfare == 1 ) then
                            call interpolation(agrid, maxG_e, 1, na, na, ia, omega)
                            V_c(i,a,j,s) = (one/(one-sigma))*(( maxC_e )**(one-sigma)) + beta(b)*theta*( omega*EVc_ge(ia) + (one-omega)*EVc_ge(ia+1) ) + beta(b)*(one-theta)*( omega*EVc_be(ia) + (one-omega)*EVc_be(ia+1) )
                        end if
                        
                    else
                    	V(i,a,j,s) = maxV_u
                    	g(i,a,j,s) = maxG_u
                    	c(i,a,j,s) = maxC_u
                    	d(i,a,j,s) = zero
                    	h(i,a,j,s) = zero
                        if ( switch_welfare == 1 ) then
                            call interpolation(agrid, maxG_u, 1, na, na, ia, omega)
                            V_c(i,a,j,s) = (one/(one-sigma))*(( maxC_u )**(one-sigma)) + beta(b)*theta*( omega*EVc_gu(ia) + (one-omega)*EVc_gu(ia+1) ) + beta(b)*(one-theta)*( omega*EVc_bu(ia) + (one-omega)*EVc_bu(ia+1) )
                        end if
                        
                	end if

            	end if

            end do

	end do
	!$omp end parallel do

	!print *, j

end do

!print *, "End Solving for Decision Rules "

    
!!Solving for Decision Rules 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!Solving for the Distribution of Households

!!!! Initialize Distribution of Households !!!!
mu = zero
mur = zero

!!!! Distribution of Households at Age 1 !!!!
do i = 1,ni
	e = e_index(i)
	n = n_index(i)
	b = b_index(i)
	if ( n==1) then
		mu(i,nd,1,1) = pi_erg(e)*phi_b(b)*phi*phi_e
	end if
	if ( n==2 ) then
		mu(i,nd,1,1) = pi_erg(e)*phi_b(b)*phi*(one-phi_e)
	end if
end do
    
!!!! Future Distribution of Households - Working Ages !!!!
do j = 1,(nj-1)
    do a = 1,na
        do i = 1,ni
            
            !!!! Good Credit Status !!!!
            s = 1
            muval = mu(i,a,j,s)
            if ( muval > precision_mu ) then
                gval = g(i,a,j,s)
				dval = d(i,a,j,s)
				hval = h(i,a,j,s)
                call interpolation(agrid, gval, 1 ,na, na ,ia, omega)
                
                do i2 = 1,ni
                    if ( hval > (one-precision) ) then
                        if ( dval > (one-precision) ) then
                            mu(i2,nd,j+1,2) = mu(i2,nd,j+1,2) + pi_i(i,i2,1)*muval
							!mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-omega)*pi_i(i,i2,1)*muval
                        elseif ( dval < precision ) then
                            mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + omega*pi_i(i,i2,1)*muval
							mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-omega)*pi_i(i,i2,1)*muval
                        else
                            mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + (one-dval)*omega*pi_i(i,i2,1)*muval
							mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-dval)*(one-omega)*pi_i(i,i2,1)*muval
							mu(i2,nd,j+1,2) = mu(i2,nd,j+1,2) + dval*pi_i(i,i2,1)*muval
							!mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + dval*(one-omega)*pi_i(i,i2,1)*muval
                        end if
                    elseif ( hval < precision ) then
                        if ( dval > (one-precision) ) then
                            mu(i2,nd,j+1,2) = mu(i2,nd,j+1,2) + pi_i(i,i2,2)*muval
							!mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-omega)*pi_i(i,i2,2)*muval
                        elseif ( dval < precision ) then
                            mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + omega*pi_i(i,i2,2)*muval
							mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-omega)*pi_i(i,i2,2)*muval
                        else
                            mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + (one-dval)*omega*pi_i(i,i2,2)*muval
							mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-dval)*(one-omega)*pi_i(i,i2,2)*muval
							mu(i2,nd,j+1,2) = mu(i2,nd,j+1,2) + dval*pi_i(i,i2,2)*muval
							!mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + dval*(one-omega)*pi_i(i,i2,2)*muval
                        end if
                    else
                        if ( dval > (one-precision) ) then
                            mu(i2,nd,j+1,2) = mu(i2,nd,j+1,2) + hval*pi_i(i,i2,1)*muval
							!mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + hval*(one-omega)*pi_i(i,i2,1)*muval
							mu(i2,nd,j+1,2) = mu(i2,nd,j+1,2) + (one-hval)*pi_i(i,i2,2)*muval
							!mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-hval)*(one-omega)*pi_i(i,i2,2)*muval
                        elseif ( dval < precision ) then
                            mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + hval*omega*pi_i(i,i2,1)*muval
							mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + hval*(one-omega)*pi_i(i,i2,1)*muval
							mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + (one-hval)*omega*pi_i(i,i2,2)*muval
                            mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-hval)*(one-omega)*pi_i(i,i2,2)*muval
                        else
                            mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + hval*(one-dval)*omega*pi_i(i,i2,1)*muval
							mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + hval*(one-dval)*(one-omega)*pi_i(i,i2,1)*muval
							mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + (one-hval)*(one-dval)*omega*pi_i(i,i2,2)*muval
							mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-hval)*(one-dval)*(one-omega)*pi_i(i,i2,2)*muval
                    
							mu(i2,nd,j+1,2) = mu(i2,nd,j+1,2) + hval*dval*pi_i(i,i2,1)*muval
							!mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + hval*dval*(one-omega)*pi_i(i,i2,1)*muval
							mu(i2,nd,j+1,2) = mu(i2,nd,j+1,2) + (one-hval)*dval*pi_i(i,i2,2)*muval
							!mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-hval)*dval*(one-omega)*pi_i(i,i2,2)*muval
                        end if
                    end if
                    
                end do
            end if
            
            !!!! Bad Credit Status !!!!
            s = 2
            muval = mu(i,a,j,s)
            if ( muval > precision_mu ) then
                gval = g(i,a,j,s)
                dval = d(i,a,j,s)
				hval = h(i,a,j,s)
                call interpolation(agrid, gval, 1, na, na, ia, omega)
                
				do i2 = 1,ni
                    if ( hval > (one-precision) ) then
						mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + theta*omega*pi_i(i,i2,1)*muval
						mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + theta*(one-omega)*pi_i(i,i2,1)*muval
                        mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + (one-theta)*omega*pi_i(i,i2,1)*muval
						mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-theta)*(one-omega)*pi_i(i,i2,1)*muval
                    elseif ( hval < precision ) then
                        mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + theta*omega*pi_i(i,i2,2)*muval
						mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + theta*(one-omega)*pi_i(i,i2,2)*muval
                        mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + (one-theta)*omega*pi_i(i,i2,2)*muval
						mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-theta)*(one-omega)*pi_i(i,i2,2)*muval
                    else
                        mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + hval*theta*omega*pi_i(i,i2,1)*muval
						mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + hval*theta*(one-omega)*pi_i(i,i2,1)*muval
						mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + (one-hval)*theta*omega*pi_i(i,i2,2)*muval
						mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-hval)*theta*(one-omega)*pi_i(i,i2,2)*muval
                    
						mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + hval*(one-theta)*omega*pi_i(i,i2,1)*muval
						mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + hval*(one-theta)*(one-omega)*pi_i(i,i2,1)*muval
						mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + (one-hval)*(one-theta)*omega*pi_i(i,i2,2)*muval
						mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-hval)*(one-theta)*(one-omega)*pi_i(i,i2,2)*muval
                    end if
                end do
            end if

        end do
    end do
end do
    
!!!! Distribution of Households - Last Working Age !!!!
j = nj
do a = 1,na
    do i = 1,ni
        do s = 1,ns
			muval = mu(i,a,j,s)
			if ( muval > precision_mu ) then
                b = b_index(i)
				gval = g(i,a,j,s)
                call interpolation(agrid, gval, 1, na, na, ia, omega)
				mur(b,ia,1) = mur(b,ia,1) + omega*muval
				mur(b,ia+1,1) = mur(b,ia+1,1) + (one-omega)*muval
            end if
            
        end do
    end do
end do

!!!! Distribution of Households - Retirement Ages !!!!
do j = 1,(nR-1)
    do a = 1,na
        do b = 1,nb
            muval = mur(b,a,j)
            if ( muval > precision_mu ) then
				gval = gr(b,a,j)
                call interpolation(agrid, gval, 1, na, na, ia, omega)
                mur(b,ia,j+1) = mur(b,ia,j+1) + omega*muval
                mur(b,ia+1,j+1) = mur(b,ia+1,j+1) + (one-omega)*muval
            end if
            
        end do
    end do
end do      
    

!! Solving for the Distribution of Households
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! End Outler Loop

!!!! Calculate Aggregate Capital and Labor !!!!
aggKp = zero; aggLp = zero
do j = 1,nj
    do i = 1,ni
        e = e_index(i); n = n_index(i)
        do a = 1,na
            do s = 1,ns
                gval = g(i,a,j,s)
                if ( gval < zero ) then
                    aggKp = aggKp + q(i,a,j)*gval*mu(i,a,j,s)
                else
                    aggKp = aggKp + gval*mu(i,a,j,s)/(one+rate_save)
                end if
                
                if ( n == 1 ) then
                    aggLp = aggLp + egrid(e)*gamma_j(j)*mu(i,a,j,s)
                end if
                
            end do
        end do
    end do
    
    if ( j <= nR ) then
        do b = 1,nb
            do a = nd,na
                aggKp = aggKp + gr(b,a,j)*mur(b,a,j)/(one+rate_save)
            end do
        end do
    end if
    
end do

!!!! Calculate Error in Outer Loop !!!!
error_K = abs(aggKp - aggK)
error_L = abs(aggLp - aggL)
error = maxval( (/error_K, error_L /) )
if ( switch_outerloop == 1 ) then
    iteration = iteration + 1
else 
    iteration = 1000
end if

!!!! Update Capital and Labor !!!!
aggK = rho_K*aggK + (one-rho_K)*aggKp
aggL = rho_L*aggL + (one-rho_L)*aggLp

if ( iteration == 1 ) then
    print '(3x, a7, 5x, a7, 5x, a9)', "Error K", "Error L", "iteration"
end if
print '(1x, ES10.4, 3x, ES10.4, 4x, I5)', error_K, error_L, iteration

end do

!! End Outer Loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Calculate and Display Results
    
!!!! Save the Value Functions to Calculate Welfare - Age 1 !!!!
if ( switch_welfare == 1 ) then
    open(unit=8, file='Welfare_Results_SS.txt')
    write(unit=8, fmt=*) "Welfare Results - SS Equilibrium"
    write(unit=8, fmt=*) "Scenario: ", switch_scenario
    do i = 1,ni
        e = e_index(i)
        n = n_index(i)
        b = b_index(i)
        write(unit=8, fmt='(I2, 3x, I1, 3x, I1, 3x, f16.13, 3x, f13.8, 3x, f13.8)') e, n, b, mu(i,nd,1,1), V(i,nd,1,1), V_c(i,nd,1,1)
    end do
end if  
    
    if ( switch_results == 1 ) then
    
        !!!! Initialize Values for Results !!!!
        aggY = zero; aggC = zero; aggI = zero; aggG = zero; aggD = zero; aggB = zero; avg_earn = zero
		Dpop = zero; Epop = zero; Upop = zero; Npop = zero; aggXi = zero
        aggC_l = zero; aggC_h = zero; aggC_debt = zero; aggC_save = zero
        Dpop_E = zero; Dpop_U = zero; Dpop_P = zero; Dpop_N = zero
        avg_q = zero; avg_limit = zero; pop_borrow = zero
        limit_all = zero; limit_patient = zero; limit_impatient = zero
        pop_all = zero; pop_patient = zero; pop_impatient = zero
        UI_exp = zero; ret_exp = zero; T_exp = zero; tax_rev = zero
        aggA = zero; agg_discharge = zero; profits=zero
        
        do a = 1,na
            do j = 1,nj
				do i = 1,ni
					e = e_index(i); n = n_index(i); b = b_index(i)
                    do s = 1,ns
                        muval = mu(i,a,j,s)
                        cval = c(i,a,j,s)
                        gval = g(i,a,j,s)
                        hval = h(i,a,j,s)
                        dval = d(i,a,j,s)
                        if ( s == 1 ) then
                            blim = abar(i,a,j)
                        else
                            blim = zero
                        end if
                        
                        if ( muval > precision_mu ) then
                            
                            !!!! Aggregate Quantities !!!!
                            aggC = aggC + cval*muval
                            aggB = aggB + dval*muval
                            T_exp = T_exp + y_bar*muval
                            if ( a < nd ) then
                                aggD = aggD + agrid(a)*muval
                                Dpop = Dpop + muval
                                agg_discharge = agg_discharge + dval*agrid(a)*muval
                                avg_limit = avg_limit + blim*muval
                                aggC_debt = aggC_debt + cval*muval
                            else
                                aggA = aggA + agrid(a)*muval
                                aggC_save = aggC_save + cval*muval
                            end if
                            
                            if ( b == 1 ) then
                                aggC_h = aggC_h + cval*muval
                            else
                                aggC_l = aggC_l + cval*muval
                            end if
                            
                            if ( s == 1 ) then
                                limit_all = limit_all + blim*muval
                                pop_all = pop_all + muval
                                if ( b == 1 ) then
                                    limit_patient = limit_patient + blim*muval
                                    pop_patient = pop_patient + muval
                                else
                                    limit_impatient = limit_impatient + blim*muval
                                    pop_impatient = pop_impatient + muval
                                end if
                            end if
                            
                            !!!! Results - Employment States !!!!
                            if ( n == 1 ) then
								Epop = Epop + muval
                                aggXi = aggXi + (1.0_rp-hval)*muval
                                aggXi = aggXi + xi*hval*muval
								avg_earn = avg_earn + wage*egrid(e)*gamma_j(j)*muval
								tax_rev = tax_rev + tau*wage*egrid(e)*gamma_j(j)*muval
								if ( a < nd ) then
									Dpop_E = Dpop_E + muval
                                end if
							else
								Upop = Upop + hval*muval
								Npop = Npop + (one-hval)*muval
								if ( a < nd ) then
									Dpop_U = Dpop_U + hval*muval
									Dpop_N = Dpop_N + (one-hval)*muval
								end if
                            end if

                            !!!! UI Expenditures !!!!
							if ( n == 2 ) then 
								if ( upsilon_r*wage*egrid(e)*gamma_j(j) >= upsilon_bar ) then
									y_u = upsilon_bar
								else
									y_u = upsilon_r*wage*egrid(e)*gamma_j(j)
								end if
								UI_exp = UI_exp + y_u*muval
                            end if
                            
                            !!!! Borrowing and Saving !!!!
                            if ( gval < zero ) then
                                qval = q(i,a,j)
                                pop_borrow = pop_borrow + muval
                                avg_q = avg_q + qval*muval
                            end if

						end if

					end do					

                end do

                !!!! Results for Retired !!!!
				if ( j <= nR ) then
					do b = 1,nb

                        !!!! Aggregate Quantities !!!!
                        cval = cr(b,a,j)
                        gval = gr(b,a,j)
						aggC = aggC + cval*mur(b,a,j)
                        aggC_save = aggC_save + cval*mur(b,a,j)
                        if ( b == 1 ) then
                            aggC_h = aggC_h + cval*mur(b,a,j)
                        else
                            aggC_l = aggC_l + cval*mur(b,a,j)
                        end if

						ret_exp = ret_exp + y_r*mur(b,a,j)
                        if ( a > nd ) then
                    		aggA = aggA + agrid(a)*mur(b,a,j)
                		end if

					end do
				end if

			end do

        end do
        aggXi = aggXi/Epop
        aggC_h = aggC_h/phi_b(1)
        aggC_l = aggC_l/phi_b(2)
        aggC_debt = aggC_debt/Dpop
        aggC_save = aggC_save/(one-Dpop)
        avg_earn = avg_earn/Epop
        avg_limit = avg_limit/Dpop
        limit_all = limit_all/pop_all
        limit_patient = limit_patient/pop_patient
        limit_impatient = limit_impatient/pop_impatient
        Dpop = Dpop/sum(mu)
        Dpop_P = (Dpop_E+Dpop_U)/(Epop+Upop)
		Dpop_E = Dpop_E/Epop
		Dpop_U = Dpop_U/Upop
		Dpop_N = Dpop_N/Npop
        
		Npop = Npop + sum(mur)
        aggI = aggKp - (one-delta)*aggK
		aggY = zgrid*((aggK)**(alpha))*(aggL**(one-alpha))
        
        avg_q = avg_q/pop_borrow

        profits = (one+rate)*aggK - aggA - aggD + agg_discharge + iota*aggD         !Note: AggKp and agg savings cancel each other out so I don't write them explicitely
        aggG = tax_rev - UI_exp - ret_exp - T_exp + profits
    
        print *, ""
        print *, ""
        print *, "----------------------------------"
        print *, " Steady State Equilibrium Results "
        print *, ""
        print *, "Agg Capital: ", aggK
        print *, "Agg Labor: ", aggL
        print *, "Rate: ", rate
        print *, "Wage: ", wage

    end if

!!!! Parameter Values !!!!
if ( switch_results == 1 ) then
    open(unit=15, file='Results_SS.txt')
    write(unit=15, fmt=*) "SS Equilibrium Results"
    write(unit=15, fmt=*) "Scenario: ", switch_scenario
    write(unit=15, fmt=*) ""
    write(unit=15, fmt=*) ""
    write(unit=15, fmt=*) "----------------------------------------------------------------------"
    write(unit=15, fmt=*) "----------------- Parameter Values - per Scenario --------------------"
    write(unit=15, fmt=*) ""
    write(unit=15, fmt=*) "TFP:       ", zgrid
    write(unit=15, fmt=*) "xi:        ", xi
    write(unit=15, fmt=*) "lambda:    ", lambda
    write(unit=15, fmt=*) "upsilon_d: ", upsilon_d
    write(unit=15, fmt=*) "upsilon_r: ", upsilon_r
end if

!!!! Aggregate Moments !!!!
if ( switch_results == 1 ) then
    write(unit=15, fmt=*) ""
    write(unit=15, fmt=*) ""
    write(unit=15, fmt=*) "----------------------------------------------------------------------"
    write(unit=15, fmt=*) "------------------------ Aggregate Moments ---------------------------"
    write(unit=15, fmt=*) ""
    write(unit=15, fmt='(24x, a4)') "Mean"
    write(unit=15, fmt='(1x, a8, 12x, f8.4)') "Output: ", aggY
    write(unit=15, fmt='(1x, a13, 7x, f8.4)') "Consumption: ", aggC
    write(unit=15, fmt='(1x, a12, 8x, f8.4)') "Investment: ", aggI
    write(unit=15, fmt='(1x, a9, 11x, f8.4)') "Gov Exp: ", aggG
    write(unit=15, fmt='(1x, a12, 8x, f8.4)') "Uns Credit: ", aggD
    write(unit=15, fmt='(1x, a14, 6x, f8.4)') "Bankruptcies: ", aggB
    write(unit=15, fmt='(1x, a16, 4x, f8.4)') "Credit Premium: ", 100.0_rp*( ( (2.0_rp-avg_q)**4.0_rp - 1.0_rp ) - ( (1.0_rp+rate)**4.0_rp - 1.0_rp ) )
    write(unit=15, fmt='(1x, a18, 2x, f8.4)') "Borrowing Limits: ", avg_limit
    write(unit=15, fmt='(1x, a14, 6x, f8.4)') "Unemployment: ", 100.0_rp*Upop/(Epop+Upop)
    write(unit=15, fmt='(1x, a15, 5x, f8.4)') "Participation: ", 100.0_rp*(Epop+Upop)/(Epop+Upop+Npop)
end if

!!!! Calibration Targets !!!!
if ( switch_results == 1 ) then
    write(unit=15, fmt=*) ""
    write(unit=15, fmt=*) ""
    write(unit=15, fmt=*) "----------------------------------------------------------------------"
	write(unit=15, fmt=*) "------------- Calibration Targets (Annualized 1980-2019) -------------"
	write(unit=15, fmt=*) ""
	write(unit=15, fmt='(a9, 3x, a5, 6x, a6, 19x, a4, 7x, a5 )') "Parameter", "Value", "Target", "Data", "Model"
	write(unit=15, fmt='(1x, a6, 3x, f7.3, 4x, a17, 8x, f7.3, 5x, f7.3 )') "beta_h", beta(1), "Capital-GDP Ratio", 224.07, (aggK/(4.0_rp*aggY))*100.0_rp
    write(unit=15, fmt='(1x, a5, 4x, f7.3, 4x, a17, 8x, f7.3, 5x, f7.3 )') "delta", delta, "Inv-Capital Ratio", 7.89, ((4.0_rp*aggI)/aggK)*100.0_rp
	write(unit=15, fmt='(1x, a6, 3x, f7.3, 4x, a14, 11x, f7.3, 5x, f7.3 )') "beta_l", beta(2), "Debt-GDP Ratio", 4.83, -(aggD/(4.0_rp*aggY))*100.0_rp
    write(unit=15, fmt='(1x, a7, 2x, f7.3, 4x, a17, 8x, f7.3, 5x, f7.3 )') "pi_beta", pi_beta, "Share of LF Borrowing", 41.37, 100.0_rp*Dpop_P
	write(unit=15, fmt='(1x, a5, 4x, f7.3, 4x, a22, 3x, f7.3, 5x, f7.3 )') "chi_b", chi_b, "Chapter 7 Bankruptcies", 0.88, 4.0_rp*aggB*100.0_rp
    write(unit=15, fmt='(1x, a4, 5x, f7.3, 4x, a14, 11x, f7.3, 5x, f7.3 )') "iota", iota_g, "Credit Premium", 11.73, 100.0_rp*( ( (2.0_rp-avg_q)**4.0_rp - 1.0_rp ) - ( (1.0_rp+rate)**4.0_rp - 1.0_rp ) )
    write(unit=15, fmt='(1x, a5, 4x, f7.3, 4x, a13, 12x, f7.3, 5x, f7.3 )') "y_bar", y_bar, "Gov-GDP Ratio", 5.38, (aggG/aggY)*100.0_rp
    write(unit=15, fmt='(1x, a5, 4x, f7.3, 4x, a18, 7x, f7.3, 5x, f7.3 )') "chi_s", chi_s, "Participation Rate", 73.73, 100.0_rp*(Epop+Upop)/(Epop+Upop+Npop)
    write(unit=15, fmt='(1x, a8, 1x, f7.3, 4x, a19, 6x, f7.3, 5x, f7.3 )') "lambda_g", lambda_g, "Mean Unemp Duration", 20.93, 13.0_rp/lambda
	write(unit=15, fmt='(1x, a4, 5x, f7.3, 4x, a17, 8x, f7.3, 5x, f7.3 )') "xi_g", xi_g, "Unemployment Rate", 6.32, 100.0_rp*Upop/(Epop+Upop)
    write(unit=15, fmt='(1x, a4, 5x, f7.3, 4x, a13, 12x, f7.3, 10x, a2 )') "xi_b", xi_b, "SD Unemp Rate", 11.12, "NA"
    write(unit=15, fmt='(1x, a8, 1x, f7.3, 4x, a15, 10x, f7.3, 10x, a2 )') "lambda_b", lambda_b, "SD_lambda/SD_xi", 1.62, "NA"
    write(unit=15, fmt='(1x, a3, 6x, f7.3, 4x, a6, 19x, f7.3, 10x, a2 )') "z_b", z_b, "SD GDP", 1.23, "NA"
    write(unit=15, fmt='(1x, a8, 1x, f7.3, 4x, a16, 9x, f7.3, 5x, f7.3)') "avg_earn", e_bar, "Average Earnings", e_bar, avg_earn
end if
    
!!!! Borrowing Limits by Household Type !!!!
if ( switch_results == 1 ) then
    write(unit=15, fmt=*) ""
    write(unit=15, fmt=*) ""
    write(unit=15, fmt=*) "----------------------------------------------------------------------"
    write(unit=15, fmt=*) "---------------------- Credit Limits by HH Type ----------------------"
    write(unit=15, fmt=*) ""
    write(unit=15, fmt='(16x, a7, 6x, a7, 3x, a9)') "All HHs", "Patient", "Impatient"
    write(unit=15, fmt='(1x, a7, 7x, f7.3, 5x, f7.3, 5x, f7.3)') "Limit: ", limit_all, limit_patient, limit_impatient
    write(unit=15, fmt='(1x, a8, 6x, f7.3, 5x, f7.3, 5x, f7.3)') "Popul.: ",pop_all, pop_patient, pop_patient
end if
    

end program BCI_SS